%% SCRIPT TO PROCESS CFRADIAL DOW DATA FOR MARSHALL FIRE TO PRODUCE VOLUMES AND SURFACE MAPS

close all
clear all


addpath('~/Dropbox/MATLAB/')
cd('~/Dropbox/MARSHALL_FIRE/')
%% Start by reading JSON file with fire perimtere data
fileName='marshall.json';
str = fileread(fileName); % dedicated for reading files as text
test=jsondecode(str);
figure(1);clf;
for aa=1:size(test.features.geometry.coordinates)
    figure(1);hold on;
    pnow=[squeeze(test.features.geometry.coordinates{aa})];
    plot(pnow(:,1),pnow(:,2),'r')
end

%%

% fname='cfrad.20211230_223205.983_ to_20211230_223217.858_DOW6low_SUR.nc'

files=dir('./deployment1/cfrad*.nc');



%% FOr the first pass through the data we need to group files into "volumes" defined as upward increments in the tilt angle
for ff=1:length(files);
    fname=strcat(files(ff).folder,'/',files(ff).name);
    azimuth=double(ncread(fname,'azimuth'));
    elevation=double(ncread(fname,'elevation'));
    evlist(ff)=mode(elevation);

end

%% Volumes should have incrementally increasing elevations and at least 4 elevation steps
figure(1);clf;
nvals=1:ff;
dev=[0 diff(evlist)];
sp1=subplot(2,1,1); cla; hold on;
plot(1:ff,evlist)
hold on;
scatter(nvals(dev<-1),evlist(dev<-1),'*k')
sp2=subplot(2,1,2);cla; hold on;
plot(nvals,dev)
hold on;
scatter(nvals(dev<-1),dev(dev<-1),'*k')

%% now find the dev<1 points as break points
bidx=find(dev<-1);
col=jet(50);
for bbb=1:(length(bidx)-1)
    if bbb==1
        segment(bbb).indices=1:bidx(bbb);
    else
        segment(bbb).indices=bidx(bbb):bidx(bbb+1)
    end

    plot(sp1,segment(bbb).indices,evlist(segment(bbb).indices))
end

%% Now loop over the number of segements, and add all the data from those files to one volume

% set up grid to house the data
xv=(-2*10^4):100:(4*10^4);
yv=(-1.6*10^4):100:(.5*10^4);
zv=[0:100:5000];
[xxx,yyy,zzz]=meshgrid(xv,yv,zv);
dbzmaster=nan(length(segment),length(yv),length(xv),length(zv));
timemaster=nan(length(segment),1);
for vv=1:length(segment)
    disp(vv./length(vv))
    fileidx=segment(vv).indices;
    allx=[];
    ally=[];
    allz=[];
    alldbz=[];
    for ff=1:length(fileidx)


        %%
        fname=strcat(files(fileidx(ff)).folder,'/',files(fileidx(ff)).name);
        ncdisp(fname)
        lats=ncread(fname,'latitude');
        lons=ncread(fname,'longitude');
        alt=ncread(fname,'altitude');
        vel=double(ncread(fname,'VEL'));
        velF=double(ncread(fname,'VEL_F'));
        dbzhc=double(ncread(fname,'DBZHC'));
        width=double(ncread(fname,'WIDTH'));

        azimuth=double(ncread(fname,'azimuth'));
        elevation=double(ncread(fname,'elevation'));
        ranger=double(ncread(fname,'range'));
        stime=[ncread(fname,'time_coverage_start')'];
        yyyymmdd=stime([1:4 6 7 9 10]);
        HHMMSS=stime([12 13 15 16 18 19]);
        timemaster(vv)=datenum(strcat(yyyymmdd,HHMMSS),'yyyymmddHHMMSS');

        re=6.37*10^6;

        clear xdist ydist dbznow

        for ev=1%:size(elevation)
            evnow=mean(elevation); % pull out the elevation for this scan.
            if evnow>0
                ae=(4/3).*(6.37*10^6); %earth radius
                z=(((ranger.^2) + (ae.^2) + (2*ranger.*ae.*(sind(evnow)))).^(1/2))-ae;
                hdist=  ae.*(asin((ranger.*cosd(evnow))./(ae+z)));%ae.*asin(rdist_HI.*(cosd(evnow))./(ae+z));
                [azmt,hdstmt]=meshgrid(azimuth,hdist); % create a mesh of azimuths and along ground distances
                [~,zdist]=meshgrid(azimuth,z);
                xdist=sind(azmt).*hdstmt; %x component of the horizontal distance
                ydist=cosd(azmt).*hdstmt;
                scalefactor=(re.*cosd(nanmean(lats))).*(pi./180);
                xdistlon=xdist./(scalefactor);
                ydistlat=ydist./(111180);

                %%
                addpath('~/Dropbox/MATLAB/SGP_LIDAR_CODES/')
                figure(1);clf;set(gcf,'pos',[100 100 1200 500],'color','w')
                sp1=subplot(1,2,1);cla;hold on;
                colormap(sp1,rbmapper(30,-10))
                pcolor(xdist,ydist,vel);shading flat;caxis([-10 30])
                cbh=colorbar;
                title(num2str(evnow),'fontsize',15)
                xlim([-1 3.5]*10^4);
                ylim([-2 1]*10^4);
                daspect([1 1 1])
                sp2=subplot(1,2,2);cla; hold on;
                %% CRITICAL NOTE HERE ABOUT NAMING CONVENTION
                dbzhc=medfilt2(width); %note I'm just keeping the dbzhc nomeclature for simplicity really thisis the spectrum width
                dbzhc(xdist<-3300)=nan;
                dbzhc(dbzhc<0)=nan;
                dbzhc(ydist>(0.2*10^4))=nan;
                dbzhc=medfilt2(dbzhc);

                %% add padded edges
                dbzhcx=imdilate(dbzhc,ones(10,10));
                dbzhc(~isinf(dbzhcx) & isnan(dbzhc))=0;


                %%
                pcolor(xdist,ydist,(dbzhc));shading flat;caxis([0 10])
                cbh=colorbar;
                colormap(sp2,jet(48))
                daspect([1 1 1])
                title(num2str(evnow),'fontsize',15)
                xlim([-1 3.5]*10^4);
                ylim([-2 1]*10^4);
                pause(.3)
                allx=[allx(:); xdist(~isnan(dbzhc))];
                ally=[ally(:); ydist(~isnan(dbzhc))];
                allz=[allz(:); zdist(~isnan(dbzhc))];
                alldbz=[alldbz(:); dbzhc(~isnan(dbzhc))];
                
            end
        end
    end


    %% visualize 3d data
    figure(188);clf;hold on;
    scatter3(allx(alldbz>1),ally(alldbz>1),allz(alldbz>1),10,alldbz(alldbz>1),'filled');
    caxis([0 5]);
    %%
    tic
    F1=TriScatteredInterp(allx(alldbz>1),ally(alldbz>1),allz(alldbz>1),alldbz(alldbz>1));
    toc
    %% create a scattered interpolant and meshed data

    tic
    dbzint=F1(xxx(:),yyy(:),zzz(:));
    toc
    dbzint=reshape(dbzint,size(xxx));
    dbzint(isnan(dbzint))=0;
    %%
    dbzmaster(vv,:,:,:)=dbzint;

end


%% OUTPUT THE DATA... NOTE: dbzmaster is really the spectrum width. I was in a rush so just renamed variables to keep code structure
stationlat=nanmean(lats);
stationlon=nanmean(lons);
stationalt=nanmean(alt);
save('marhsall_plume_d1_SW.mat','dbzmaster','timemaster','xxx','yyy','zzz','stationlat','stationlon','stationalt','-v7.3')

